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Abstract 

In this article, we study bivariate polynomial interpolation on the node points of degenerate 
Lissajous figures. These node points form Chebyshev lattices of rank 1 and are general¬ 
izations of the well-known Padua points. We show that these node points allow unique 
interpolation in appropriately defined spaces of polynomials and give explicit formulas for 
the Lagrange basis polynomials. Further, we prove mean and uniform convergence of the 
interpolating schemes. For the uniform convergence the growth of the Lebesgue constant 
has to be taken into consideration. It turns out that this growth is of logarithmic nature. 
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1. Introduction 

A by now well-established point set for Lagrange interpolation on [—1, l] 2 is given by the 
Padua points [9]. This set of points allows unique interpolation in the space II 2 of bivariate 
polynomials of degree n and has a few outstanding properties: it can be characterized as a set 
of node points of a particular Lissajous curve [3], as an affine variety of a polynomial ideal J6] 
and as a particular Chebyshev lattice of rank ira- Moreover, the Lagrange interpolant can 
be computed in a fast and efficient way and the asymptotic of the corresponding Lebesgue 
constant is of log-squared type mm- 

Using affine mappings of the square [—1, l] 2 , the Padua points can also be used for 
interpolation on general rectangular domains m- However, for highly anisotropic rectangles 
the particular structure of the Padua points is not so well adapted. In this case it is more 
favorable to use interpolation nodes that reflect different resolutions along the axes of the 
rectangle. To obtain more flexibility for anisotropic domains, it is therefore reasonable to 
consider generalizations of the Padua points. 
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In [T21126] , a framework for approximation and interpolation on general Chebyshev lat¬ 
tices was developed. This framework contains multidimensional anisotropic lattices and the 
Padua points are included as a special case. But, in this framework the interpolation of 
functions is in general only available in an approximative way. This approximate interpola¬ 
tion is known as hyperinterpolation [28j and is frequently used in applications. For trivariate 
Lissajous curves, it is studied in the recent work [5]. 

A different way of generalizing the theory of the Padua points was recently presented 
in da In this article, the node points of non-degenerate Lissajous curves were used as 
interpolation points. These nodes turned out to be special Chebyshev lattices of rank 1. 
In particular, in this theory also anisotropic point sets can be chosen. However, since the 
generating curve of the Padua points is a degenerate Lissajous curve, the Padua points are 
not directly included in the theory developed in ra. 

The goal of this article is to develop an interpolation theory for node points of degenerate 
Lissajous curves that contains the Padua points as a special case. To this end, we rebuild 
the interpolation theory developed in da for degenerate Lissajous curves. While most of the 
results carry over, some technical aspects in the proofs differ considerably. This is mainly 
due to the different geometric properties of the underlying curves and its interpolation nodes. 

The search for favorable node points in multivariate polynomial interpolation has a long¬ 
standing history. We refer to the survey articles [T8l ,19] for a general overview. Beside the 
Padua points, the points introduced by Morrow and Patterson [25] and Xu Ed (see also 
[20]) are of special importance for the theory presented in this article. 

We start our research by studying the node points LD np of degenerate Lissajous curves. 
It turns out that these node points can be characterized in several ways and, in particular, 
as Chebyshev lattices of rank 1. As in the non-degenerate case, the node points LD„ )P 
can be used as quadrature rules on [—1, l] 2 for integrals with a product Chebyshev weight. 
Compared to non-degenerate Lissajous curves, the node sets LD np in the degenerate setting 
are smaller, contain more asymmetries and include two vertices of the square [—1, l] 2 . 

The main results of this article can be found in Section 4 and 5. In Theorem [7] we 
prove that LD njP allows unique polynomial interpolation in a properly defined space n 2, p of 
bivariate polynomials. We will derive an explicit formula for the corresponding fundamental 
Lagrange polynomials. This explicit formula is very similar to the one known for the Padua 
points and allows to compute the interpolating polynomial in a simple and efficient way. 

Whereas in [15] stability and convergence was studied only numerically, in this article 
we investigate the convergence of the interpolation scheme also in an analytic way. For 
continuous functions / we show mean convergence of the Lagrange interpolant in the L r - 
norm. For the Xu and the Padua points it is known that the Lebesgue constants grow as 
0( In 2 n) (cf. [HE] [13]). We will confirm a similar log-squared behavior also for the Lebesgue 
constant of the general interpolation scheme considered in this article. We conclude this 
article with some numerical experiments that confirm the convergence results and illuminate 
the role of the parameter p in anisotropic setups. 
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2. The node points of degenerate Lissajous curves 

In this work, we consider Lissajous figures of the type 

7 n,p ■ [—1, l] 2 , In,pit) = (cos(nt), cos((?r + p)t) S j, (1) 

with positive integers n and p such that n and n + p are relatively prime. The curve is 
27T-periodic but doubly traversed as t varies from 0 to 27 r. For this reason, y niP is referred 
to as degenerate Lissajous curve and we can restrict the parametrization of the curve to the 
interval [0, 7r], The points 7n, P (0) = (1, 1) and 7 njP (7r) = ((— l) n , (— l) n+p ) denote the starting 
and the end point of the curve 77 . p , respectively. A classical reference for the characterization 
of two-dimensional Lissajous curves and its singularities is the dissertation [7J of Braun. In 
recent years, Lissajous curves are particularly studied in terms of knot theory 121 EH- 
If we sample the curve 7 n>p along the n{n + p) + 1 equidistant points 


tk ■ — 


7T k 

77 ,( 77 , + p) ’ 


k 


0 ,..., n(n + p), 


in the interval [ 0 , 7 r], we get the following set of node points: 

LD n ,p ■ ^7fi,p(^fc) • t 0,... , 77 .( 77 . T p) j' • (2) 

To derive particular properties of the Lissajous curve 7 „ )P and the set LD np it is easier to 
characterize the curve 7 ^ with help of the Chebyshev polynomials T n {x) = cos (77 arccosx) 
of the first kind. Based on this relation, the following properties of 7 n)P are proven in 
[17] [Section 3.9] and [2T] . We use the notation 


zj? := cos 



77 G N, k — 0, ..., n, 


to abbreviate the Chebyshev-Gaufi-Lobatto points. 


Proposition 1. If n and n + p are relatively prime, the Lissajous curve ^ n - p , t € [0,7 r] ; 
corresponds to the (plane) algebraic curve 

C n , P ■= {(x, y) e [-1, l] 2 : T n+p (x ) - T n (y) = 0}. (3) 

The curve 7 n )P {t), t G [ 0 , 7 r], has d i +p-L(n-i) or ^ nar y self-intersection points in the interior 
of the square [—1, l] 2 . They are given as 

{ i = l,...,n + p- l | 

k" +p ,y): J — 1,.... n — 1 1 (4) 

i + j = 0 mod 2 J 


and can be arranged in two rectangular grids. 
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Figure 1: Illustration of the degenerate Lissajous curve 73 , 2 ! its node points LD 3 2 and the corresponding 
index set Tg 2 according to the characterization (111 in Proposition [ 2 ] 


Using ([lj or (|3]), it is easy to see that the Lissajous curve 7 „ iP touches the boundary of 
the square [—1, l] 2 at exactly 2 n + p points. We collect these boundary points in the set 


(5) 



i — Q,...,n + p ) 

1 I 

f ie{0,n+p} ) 


: j G {()•//} ' 

U 

(zt +p ,z^ : j = 1 ,... ,n — 1 \ 

r ; 

i + j — 0 mod 2 J 

1 1 

{ i + j = 0 mod 2 J 


We can further divide the boundary points LD°p in vertex and edge points: 

LD ”? := {7nj>(0),7nj>M}> 

LD* := LD™p \ LD™". 

Now, we get the following characterizations for the node set LD np : 

Proposition 2. If n and n + p are relatively prime, the set LD n)P contains U+p+LU+L 
distinct points and is the union of self-intersection and boundary points of the curve jn !P , 

i.e. 

% ' = - 0, . . . , Ti p 

LD. W = LD“ ULD“= ( h.” +P .y) • 3 = 0 , ■■■’,» }. (6) 

i + j — 0 mod 2 
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LD„ p can be arranged in two rectangular grids 


LD” := ' (*” +, ',4) : J=0,...,n 


U lP 


LD:, p := { P”+V" 


n+p n 
i J 


z = 0,..., n + p 


i,j even 

i = 0,... ,n + p 
j = 0,.. ■, ra 
zj odd 


Further, introducing the index sets 


r n,p ■= \ ihj) e N 2 


n + p n 


+ -<![>, 


r^:=r BiP u{(0,n)}, 


the node set LD niP can be characterized as 


L D n _ n 


iP { (+n++n+p)’ ^m-h?(n+p)) ■ J) ^ 


(7) 

( 8 ) 


(9) 

( 10 ) 


( 11 ) 


PROOF. Clearly, { (zl + ^ (n+p)1 z^ n+Kn+p) ) : (i,j) G T^} is a subset of LD„ jP . For (z,j) G 

r^ p , the integers k — in + j(n + p) and k' = |zn — j(n + p) are nonnegative and less or 
equal to n(n + p). Evaluating the Lissajous curve 7 U)V at t k and t k ', we observe that the 
points 7 n ,p(tk) and 7 n ,p(tk') are equal. Further, since n and n + p are relatively prime, k and 
k' coincide if and only if i = 0 or j = 0. Therefore, for all positive integers i , j G I with 
< 1, we obtain a pair of distinct points t k ,t k ' G [0,7r] such that 7 n ,p{tk) — 7n,p+')- 

The total number of distinct pairs is given by , Therefore, by Proposition [lj 

these pairs describe exactly all self-intersection points of the curve 7 n , p , i.e. 


LD 


int 

n,p 


I n+p n \ 

\ Z in+j(n+p)i Z in+j{n+p) J 


+ 3 ) e N 2 , 


i j 

- 1 — < 

n + p n 



Further, the numbers k = in + j(n + p) with (i,j) G T^ p and i = 0 or j = 0 correspond 
precisely with the 2 n + p boundary points of 7 ^, i.e. 


LD 


out 

n,p 


n+P n 

^in+j {n+p) 1 ^in+j(n+p) 


) : (hi) e r n,p) i = 0 or 3 = 0 


Therefore, 


LD 


int 
n ,p 


U LD 


out 

n,p 


n+p n 

*in+j(n+p)i ^in+j(n+p ) 


) : (i, j) 6 !+} C LD 


y n,p • 


Since 2|LD“p| + |LD° U p| = n(n + p) + 1, we even have equality in the last formula. This 
together with the definitions of LDPl and LD°“p in Q and ([5]) implies equation ([ 6 ]) as well 
as equation (0- Finally, (J7|) and (J8| follow from (jbj). □ 
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Regarding the cardinality of the subgrids LD^ p and LD™ p , we can distinguish between 
three cases depending on whether n and p are even or odd integers. The same holds for 
the location of the second vertex 7n !P (7r) = ((—l) n , (—l) n+p ), whereas the first vertex is 
always given as 7„ iP (0) = (1,1). Using the formulas (JT]) and ([8]), we compute the different 
cardinalities and list them in Table [TJ The respective cases are also illustrated in Figure [2] 


Table 1: Number of points in the different LD sets. 


|LD njP 

(n+p+l)(n+l) it T-vint 

2 1 * J ^'n,p\ 

(n+p—l)(n—1) 

2 

|LD°^| 


= 2 n + p 


Case (a): 

n even, p odd 

Case (b): n odd, p odd 

Case (c): 

n odd, p even 

|LD^| = 
|LD- p | = 

n+2 n+p+1 

2 2 
n n+p+1 

2 2 

(1,-1) 

It p\b 1 _ n+1 n+p+2 

It T^w I _ n+1 n+p 

L1J n,p 2 2 

7 n,p( 7r ) — ( — 1) 1) 

|LD^ p | = 

|ld; p | = 

7 n,p( 7r ) — 

n+1n+p+1 

2 2 
n+1n+p+1 

2 2 

+ 1,-1) 



Figure 2: Illustration of degenerate Lissajous curves 7 njP and its node points LD„ iP according to the cases 
in Table [lj The points in the subgrids LD^ p and LD^ p are colored in blue and white, respectively. 

To simplify the notation of the integers k in ([2]) that describe the same point A G LD niP , 
we introduce on Z the equivalence relation 1 rs./ P by 

k L ~’ P k' & 7n,p(4) = 7 n,p{tk')- (12) 

In this way, we obtain for each A G LD njP a unique equivalence class [.A] and we say that 
k G [A] if 7 n ,p(tk) = A. By the argumentation in the proof of Proposition [2j there is exactly 
one 0 < k < n(n + p ) in the equivalence class [A] if A G LD°“p and exactly two if A G LD“ p 
is a self-intersection point. 
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Remark 1. In this article, the parameter choice p — 1 corresponds to the well-known Padua 
points studied extensively in PI H El El El EDI E]. For the Padua points, one can differ 
between four families of interpolation points by considering ninety degree rotations of the 
set LD nj i. Similarly, one gets four families of interpolation nodes for general p. In the cases 
(a) and (b) considered in Tableland Figure [2] this is also done by rotating LD njP by ninety 
degrees. In the case (c) one has to combine ninety degree rotations and reflections with 
respect to one of the coordinate axes to obtain the four families. The respective generating 
Lissajous figure is a rotated and reflected version of 7 niP . The interpolation theory developed 
in this paper can be applied to all four families of nodes. For simplicity, we will only consider 
the family LD niP generated by r y ntP . 

Remark 2. The set LD np is a two-dimensional Chebyshev lattice of rank 1. According 
to the notation given in [T2\ EE], the respective parameters of the Chebyshev lattice are 
d = n{n + p) and z = [n, n + p\. The following sections will show that lattices of this form 
allow unique polynomial interpolation. 


3. Lissajous node points and quadrature 

The point set LD np can be used for interpolation purposes as well as for quadrature rules 
on [—1, l] 2 . In this section, we study first a quadrature formula based on point evaluations 
at the set LD np . Of essential importance in our considerations is the operator 

E 7 :C([-1,1] 2 )^C([0,7t]), E 7 /(t) =/( 7niP (t)), t G [0,7r], 

that restricts a continuous function on [— 1 , l ] 2 to the Lissajous trajectory 7 n>p . It is clear 
that E 7 maps bivariate algebraic polynomials to even trigonometric polynomials on the 
interval [ 0 , 7 r]. 

To specify the spaces of bivariate polynomials we introduce 

n 2 = span {Ti{x)Tj{y) : i + j < n}, (13) 


where T t (x) = cos(f arccosa;), as before, denote the Chebyshev polynomials of the first kind. 
It is well-known (cf. |3l| ) that {Tj(x)Tj(y) : i + j < n} forms an orthogonal basis for the 
space II 2 with respect to the inner product 


if, 9) 


7T 


I t J J{x,y)g(x,y)-j= 


X‘ 




r 


idicdy. 


(14) 


With the normalization 

Ti(x) 


1 , if i — 0 , 

V2Ti(x), if i ^ 0 , 


we obtain the orthonormal basis {Ti(x)fj(y) : i + j < n} of II 2 . 

The first auxiliary result shows that for a large class of bivariate polynomials P the 
restriction E 7 P can be used to convert a double integral into a one dimensional integral of 
a trigonometric polynomial. 
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Lemma 3. For all bivariate polynomials P with (P,Tk( n + p )(x)Tk n (y)) = 0, 6 N, the 

following formula holds: 


»1 rl 


7T 


P{x,y) 


'-iJ-i 


\/l - \/l - y 2 


1 1 r 

cLrch/ = — / E 7 P(t)df. 


7T 


(15) 


PROOF. Lennna[3]is a slight generalization of [3] Lemma 1]. With minor changes, the proof 
follows the same lines of argumentation as the proof of [3], Lemma 1] and [15] Lemma 1], □ 

Next, we consider the spaces of bivariate polynomials corresponding to the index sets 
r n , p and T r np introduced in 


and (10): 


n 


n,p 


:= span {T i (x)T j (y) : (■ i,j ) e T njP }, 


n S := span {Ti{x)Tj{y) : (i,j) e T^ p }. 

By the characterization of the set LD np given in ( fTTj ) it follows immediately that 

(■n + p+ l)(n + 1) 


(16) 


dim = |r^ p | = |LD n>p | = 


dim n^ p = |r n , p | = \r^ p \ -1 = 


(n + p+ 1 )(n + 1) 


- 1. 


To characterize the range of the operator E 7 , we need particular spaces of trigonometric 
polynomials on [0,7r]: 


N 


:= < q(t) = ^2 a m COS (mt) : a m G M, t G [0, tt\ > . 


m =0 


Lemma 4. The operator E 7 defines an isometry from the polynomial space Iff p (n'f p ) 
equipped with the inner product (-, ■) given in (14) into the trigonometric space II®/ + % l 

i r 

(II®( n+p ), respectively) equipped with the inner product ( 51 , 92 ) — 


PROOF. An orthonormal basis of fl^ is given by < Ti(x)Tj(y ) : (i, j) G r„ iP >. The image 


e i,j(t) ■= e 7 pi(x)Tj(y))(t), (i,j) g r n , p , 

of this orthonormal basis is explicitly given by 


= < 


( 1, 

y/2 cos (int ), 
v^cos (J{n+p)t ), 

2 cos (int) cos (j(n + p)t ), 
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if (i,j) = (0,0), 

if j — 0, i < n + p, 

if i — 0, j < n, 

for all other (i,j) € r n , p . 


( 17 ) 











From the definition ([9]) of the index set we obtain in+j(n+p) < n(n+p ) for (i,j) G r n p . 
Therefore, all e^j are trigonometric polynomials in the space ^d e n ( n +p )_1 an d we can conclude 
that E 7 maps Il 2 p into the space ff®( n+p) _ 1 . 

Since r^ p = T niP U {(0, n)} and eo, n (t) = E 7 (T n (y))(t) = v^cos (n(n + p)t) is a trigono¬ 
metric polynomial of degree n(n + p) , we can also conclude that E 7 maps L[ 2 ; p into the space 

TT C 

n(n+p) 

For Pi,P 2 ^ n^’p, the product PiP 2 is a polynomial in the space n 2 n+2p _ 2 and satisfies 
(P 1 P 2 , T n+p (x)T n (y)) = 0. Therefore, by Lemma [ 3 J the set {e^ : (i,j) G T^ p } is an or¬ 
thonormal system in L1C +p ^, and thus, E 7 an isometry from II 2 ; p into the space n®, n+p x. By 
the same argumentation, the set {e t j : (i,j) G r n>p } is an orthonormal system in 
This proves the isometry also for the space II 2 . □ 

For points A G LD niP , we introduce the quadrature weights 


w A : = 


and obtain the following quadrature rule for the node set LD n p : 

Theorem 5. For all polynomials P G n| ra2 p the quadrature formula 

1 1 
\/l - x 2 y/l - y 2 

is exact. 


dxdy = ^ w a P(A) (18) 

w4.GLD n) p 



2 n(n + p) ’ 

77 ,( 77 . + p) ’ 

2 

77(77 + p) ’ 


HAeLD^, 
HAe LD^J*, 
if .4 6 LD“, 


PROOF. For all even trigonometric polynomials q G n 2n ( n +p)-i the following composite 
trapezoidal quadrature rule is exact (see [ Ml Chapter X]): 

1 r 1 /l n(n+p)-l \ 

? I q(t)dt = (2 ?(0) + x 9 (4) + 2« W J ' 


Since by the definitions (13) and (16) of the polynomial spaces n 2n+2p _ 


n 2 

ll 2n,2p 


r TT 2 

C- il 2 n+2p-l 


an d n ln,2p 


T, 


n+p(x)T n (y) , Lemma [3] yields the identity 


and II 2 n 2 p we have 


1 r 1 r 1 11 1 r 

—a / P(x, y) , . dxdy = — / E 7 P(t)d£ 

vrTrio 7 w 


For a polynomial P G fl 2 n 2 , the image E 7 P is by a similar argumentation as in the proof 
of Lemma |4j an element of n 2 n (n+ P )-v Therefore, combining the two identities above and 
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using the definition ([2]) of the points LD np as well as its characterizations in Proposition [2j 
we get the quadrature formula 

— / / P(x,y)—. —p dxdy = V waP(A). 

□ 


Remark 3. For the case p = 1 of the Padua points, Theorem [5] is proven in a similar 
way in [3]. Quadrature formulas similar to (18) exist also for other related point sets as 
the Xu points (see 12011251 EH). The operator E 7 was already introduced in [IS] to describe 
the relation between algebraic polynomials on [—1, l] 2 and trigonometric polynomials on the 
Lissajous trajectory. However, due to the geometric differences between degenerate and non¬ 
degenerate Lissajous curves, the range of the operator E 7 consists of different trigonometric 
polynomials in the two cases. 


4. Lissajous node points and interpolation 

In this section, the object of study is the following interpolation problem: for given data 
values G i at the node points A G LD np , find a unique interpolating polynomial C U)P f 
in [—1, l] 2 such that 

C n .pf{A) = /a holds for all A G LD niP . (19) 

In the bivariate setting it is a priori not clear which polynomial space has to be chosen 
for the interpolation problem (19). Since dimn 2, jj = |T^J = |LD njP |, our primary choice 


n,p l 


IS 


the polynomial space H^p defined in (16). To prove the uniqueness of (19) in n 2, (j, we 


study first an isomorphism between H 2, p and the subspace 

n nfn+p) '■= { 9 e n n(n+p) : ?(*fc) = QM for all fc, k' with k 
of the even trigonometric polynomials n®^ n+ 


n,p i 


LD n?p 

rv_/ 


k' 


( 20 ) 


Theorem 6 

n c,L 


The operator E 7 is an isometric isomorphism from, n 2 ’); onto the subspace 


n(n+p) ’ e< l u wped with the inner product given in Lemma f 


PROOF. By Lemma 

'~fn,p(fk 


e LD£ is a se 


we know that E 7 is an isometry from n 2, p into the space n®,„, If 


±o cxn 10V11U.U 1J iiurn ±± n p 11WIK ±1 n ( n+p ) • 

f-intersection point of the Lissajous curve 7 niP , then 7 n ,p(tk') = ln,pifk) 
holds for a 0 < k' < n{n + p), k' ^ k and the values E 7 P(ffc) and E 7 P(t*./) coincide. This 
property is exactly encoded in the definition (20) of the space H c n ^ i+p y This implies that 
the operator E 7 maps n 2 ’jj into the subspace 


n(n+p )' 


Now, it suffices to show that the dimensions of n 2 ’ p and Ilh^, coincide. Then, we 
immediately obtain the surjectivity and, hence, the bijectivity of E 7 . To this end, we consider 


in n e n ( n+p) the Dirichlet kernel (see 

n(n-\-p) — 1 


X Section 3]) 


D n ( n +p){t) := 1 + 2 ^2 cos (kt) + cos (n(n + p)t) = 


sin [n(n + p)t) cos | 


k =1 


sin \ 
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The translates n ^ +p) -P n (n+p) {t — 4), k = 1,... ,2n(n + p), form a basis for the space of 

trigonometric polynomials of the form q(t) = Ylm=t) P> <lm cos ( m t) + 6 m sin(mt). 

Therefore, for the space n®^ n+p N of even trigonometric polynomials we get 


D 


k 

n(n+p) 


(t) ■ = 


D n (n-\-p) (t) 

n(n+p) ’ 

^n(n+p) ^fc) 

n(n+p) 

Dn(n-\-p) 7I ") 

n(n+p) ’ 


if k — 0, 

if k — 1,..., n(n + p) — 1, 
if k = n(n + p), 


as a fundamental basis of Lagrange polynomials with respect to the equidistant points 4, 
k — 0,, n(n + p), i.e. 


D t(n+ P ) (4') = 0 <k,k'< n{n + p). 

i.L 

i(n+p) ’ 


Since not all D^ n+ , are contained in the subspace II^ +p ^ we introduce for A G LD np the 
linear combinations 


UW := £ 


Dnin+p)W' 


( 21 ) 


0 <k<n(n-\-p)\ 
ke[A] 


where [*4J denotes the equivalence class introduced in (12). The trigonometric polynomials 
l A are elements of the space fI^ + > and 


4i(4) 


1 if k e [A\, 

0 if k ^ [A]. 


( 22 ) 


Since DK + s, 0 < k < n(n + p), form a basis in the space IF , the system {l A : 
A G LD np } is a basis in the linear subspace II 

dim ( n n(t+J = I LD «-pI = dim ( n n,p)- 


e,L 
n(n+p) * 


This implies the desired equality 

□ 


Finally, we prove the uniqueness of the interpolation problem (19) in IT^ p and give explicit 
formulas for the corresponding fundamental Lagrange polynomials. As a final auxiliary tool, 
we need the reproducing kernel K T np : R 2 xR 2 —y M of the polynomial space 11^ p . It is defined 
as 

K n, P {x A ,y A ]x B ,y B ) := T i ( x A )f i (x B )f j (y A )f j (y B ). (23) 

(*j)er Lp 

With K n p we denote the respective reproducing kernel for the subspace 11^ . 


Theorem 7. For A = (. x A ,y A ) G LD np; the polynomials L A := E 1 l A are given as 


La{x,v) = ™A yK^ P (x,y,x A ,y A ) ~ \T n (y)T n (y A ) 
= U! A [K n> p(x,y,x A ,y A ) + \ f n {y)T n (y A ) 


(24) 
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and form the fundamental Lagrange polynomials in the space B‘f r p with respect to the point 


set LD np? i.e. 


La{&) — &a,Bi A, B G LD n]P . 

The interpolation problem (19) has a unique solution in B‘f T p given by 

£n, P f(x,y)= Y fALA(x,y). 

*4ELD n)P 

PROOF. Theorem 6jand property (22) of the basis functions l a G H c n ^ n+P ) introduced in (21) 
imply that the polynomials La = E^ 1 1a satisfy La{B) = 5a,b for B G LD np . Moreover, 


since the system {l a ■ A G LD rap } is a basis in LT y ', r f n+p ^, Theorem |b 
system {La ■ A G LD„ p } forms a basis of Lagrange polynomials in t 

Finally, we compute the decomposition of the trigonometric polynomials l a in the or- 


implies that also the 
lie space n^’ p . 


thonormal basis e^j. Then, using the inverse E“ , we get the representation (24) for the 
Lagrange polynomials La- 

For the calculations, we use the characterization ( fTT| for points A G LD rap , i.e. A = 
(x A , Pa) = ( Z k +P i z k ) with k = rn + s(n + p) and (r, s) G T^ p . 

We take first a look at the self-intersection points A G LD“ p . From the proof of Propo¬ 


sition [2j we know that the second integer 0 < k! < n(n+p) representing the self-intersection 
at A is g 

l A® = 


point A is given by k! = |rn — s(n+p) |. In this way, for the functions l a, A G LD™ p we get 


rtPTp) + D n(n +r )W 


n(n + p) 
2 

n{n + p) 


n(n+p) 

1 + E (2 - n(n+p )) 
m= 1 

n(n+p) 


COS 


krrnv 
n(n + p) 


+ cos 


k'mn 
n{n + p) 


cos mt 


x -> rmn smir 

1 + 2_^ ( 2 ~ 5m, n ( n+p )) COS---COS-COS mt 


m= 1 


n + p 


n 


where 5 ipj denotes the Kronecker delta, i.e. 5 tpj — 1 if i — j and 5 u j = 0 otherwise. Analogous 
calculations for the vertex points A G LD^ and the edge points A G LD® d | e yield the 


n,p 


same 


formula for Ia differing only in the weight wa- So, for general A G LD njP , we get 


n(n+p) 

Ia(B) IXA J 1 T ^ ^ (2 <^m,n(n+p)) 

m= 1 


rm'K smTT 
cos-cos-cos mt 


n + p 


n 


Now, using formula (pT|) for the basis functions e iv we get after a small calculation the 
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following formula for the coefficients ( l A ,ei,j) = \ Jq ^A(t) e i,j(t)dt: 


U-U, 

WA V2 COS ^ 7r ’ 
(U,eij) = { w a V2cos^, 
™aV2 cos|g, 
w A 2 cos cos —, 

n+p n 7 


if (b j) = (0,0), 

if i = o, j = n , 
if i = 0, j < n, 
if i ~f~ 0,3 = 0, 

for all other (i,j) G T % , 


= wa 


\T n (yj i), 


if i = o ,j = 71, 

Ti(x A )Tj(y A ), if (i,j) G T„ iP . 

Therefore, the decomposition of l A in the basis e h j can be written as 

U{t) = -Y^n{.y A )eo,nit) + Y w ATi(x A )T j (y A )e iJ (t). 

(iyj)€Fn,p 


Now, the inverse mapping E 1 and the definition (|23l) of yield 


n ,p 


L A (x, y) = E 7 l A (x, y) — w A ( K np (x, y; x A , y A ) - -T n (y)T n (y A ) 


□ 


1 


0.5 


0.8 

0.6 
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-0.5 


-1 - 
-1 


-0.5 


0.4 
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0.5 
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0.2 
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- 0.2 



(a) Lagrange polynomial L A , A = (z^jZq). 


(b) Lagrange polynomial L A , A = (z]q , zf). 


Figure 3: Space localization of the Lagrange polynomials L A £ IIg ’2 for different A £ LDg .2 


Formula (24) in Theorem [7] allows to compute the fundamental Lagrange polynomials L A 
in an explicit way. Two examples of L A are illustrated in Figure |3j Moreover, using formula 
(24) we can rewrite the interpolating polynomial £ ntP f(x,y) in terms of the orthonormal 

In this way, we obtain the representation 


Chebyshev basis {7',(./•) /)(//) 


(hJ) e r n, P J 


£n, P f(x,y)= Y c i,jTi(x)fj(y ) 

(*j)er Lp 
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with the Fourier-Lagrange coefficients = (C np f, T t {x)T 3 (y)) given by 


w A f A T i (x A )T j (y A ), if (i,j) G r n)P , 
»4.GLDn,p 

^ w AfAT n (yA), if (i,j) = (0,n). 

*4ELD n , p 


( 25 ) 


The computation of the coefficients Cjj in (25) can be formulated more compactly using 
the matrix notation 


C n ,p = {T x (LB n>p ) D / (LD n>p ) T y (LD n , p ) r ) © M n , p , (26) 

where © denotes pointwise multiplication of the single matrix entries. Here, the matrix 
C Hj p = ( Cij ) G M( n +p) x ( n+1 ) contains the coefficients 


C i,j ~ 


= f (CnJ,fi(x)Tj(y)), if (i,j) e rjj,, 

0, otherwise. 


The given data f A and the weights w A are arranged in the diagonal matrix 
D/(LD n>p ) = diag (w A f A , A G LD n , p ) G M I LD —' x l LD —I. 
Further, all evaluations Ti(x A ) and Tj(y A ) are collected in the matrices 


T x (LD np ) — 


/ 

Tq(x a ) 

\ 

j T y (LD n , p ) — 

/ 

T 0 (y A ) 

\ 

V 

Tn+p— 1 (*©l) 



V 

T n {yX) 

/ 


V 

«AGLD n? p 




«AGLDn,p 



Finally, the mask M.„ iP = {rn l)3 ) G M( n+ P) x ( n+ b i s given by 

( i, if (M')Gr BiP , 

V 2 ; if (i,j) = (0,n), 

l o, if (ij) £ r£p. 

For H C M 2 , also the point evaluations C ntP f(B ) of the interpolation polynomial can be 
written compactly as vector-matrix-vector product 


C n , P f(B) = T x (B) T C n , p T y (B) 


(27) 


with T X (B) T = (f 0 (x B ), • • • , T n+p _i(x B )) and T y (B) T = (f 0 (y B ), ■ ■ • , T n (j/ B )^. 

Remark 4. The formulation of Theorem[7]is almost identical to the formulation of Theorem 
3 in p2]- The main structural difference between the two results lies in the form of the index 
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set r£ p and in the weights w A . In the case of non-degenerate Lissajous curves, the index 
set T% p is asymptotically 4 times larger than for degenerate Lissajous curves. Also, in the 
non-degenerate case there are no vertex points. Structural differences can be also found in 
the technical aspects of the respective proofs. This is again due to the fact that the operator 
E 7 maps onto different spaces of trigonometric polynomials in the two cases. 

For the Padua points (p = 1), Theorem [7] was proven in a different way. The uniqueness 
of the interpolation problem was derived using an ideal theoretic approach (see 0 ), whereas 


formula (24) was proven using an explicit formula for the reproducing kernel. 


Also, the matrix formulations in (26) and (27) are very similar to the formulations in 
. Here, the main difference lies in the form of the mask M„.„ which encodes the index 


set T 


L n,p 

. For the Padua points the mask M np is an upper left triangular matrix (cf. [llj j 


Remark 5. If we substitute the index set in (10) by 


n,p 


r n,p : = r„,pU{(n+p,0)} 

all the results of this paper can be proven in an analogous way also for the altered index set 
rj p . In this way, we obtain in a respectively altered version of Theorem [t] the polynomials 


L A (x, y) = w A ( K n j,(x,y;x A ,yj i) + -T n+p (x)T n+p (x A ) 


A G LD 


n,p, 


as fundamental Lagrange polynomials in the space n 2 ’jj := span {Ti(x)Tj(y) : (■ i,j ) G T^ }. 
In particular, we can see that it is not possible to speak of a "natural" polynomial space for 
interpolation on the point set LD niP . For the interpolation problem (19), several reasonable 
choices are possible. 


5. Convergence results for the interpolation polynomials 

In this section, we study the convergence behavior of the interpolating polynomial C HiP f 
to a given function / G C([— 1, l] 2 ) if n gets large. In particular, for 1 < r < oo we prove 
mean convergence of the Lagrange interpolation in the L r -norms 

Further, we will give an upper bound for the growth of the Lebesgue constant. 

The proof of the mean convergence as well as the estimate for the Lebesgue constant 
depend on the following forward quadrature sum estimate. In the upcoming results, we want 
to keep the parameter p as variable as possible. For this reason, we will take particular care 
that the constants in the estimates are independent of p. 

Lemma 8. For 1 < r < oo, the inequality 


£ 

-4eLD n 


w A \P(A)\ r <C r - 


7T 


LJ.A^vt 


X‘ 


y/i 


y 


-Ax&y 


holds for all polynomials P G A 2 n ' p with the constant C r = (r + l) 2 e 2 (l + A-)" 
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( 28 ) 













PROOF. Based on an idea given in |3T| , we use a univariate inequality |22, Theorem 2] 
to estimate the quadrature sum. Applied to an even trigonometric polynomial q(t) = 
ELo Ck cos this inequality has the form 

£ l«(»i)| r < (n + £ mritt (29) 

with 0 < 9 0 < 9i < ■ ■ ■ < 9 m < 27t and £ = min(#i — 9 0 ,... ,9 m — 9 m _ i, 2n — 9 rn + 9 0 ) > 0. 
To estimate the quadrature sum, we use the grid characterization of the set LD„ iP given in 
(Jr]) and (J8]) . To simplify the calculations, we will only consider the case where n is odd and 
p is even, i.e. case (c) in Table [T} The estimates for the cases (a) and (b) in Table [I] follow 
analogously. We obtain 


£ wa\p(a)\ 


f n-\-p — 1 n — 1 
2 2 


.AeLD„ 


2 n(n + p ) 


£ £(2-M(2-^o)l^(*-r.4y 

i =0 j =0 


n+p— 1 n — 1 
2 2 


+ £ £(2-fe + i,„ +P )(2-fc J+1 ,„)|P(4y,.4 + i)l r ] • 

2=0 j =0 


Now, for every fixed 0<i< n + p — 1, we adopt inequality (29) to the even trigonometric 
polynomials P(z 2 + p ■, cost) and P(^ 2 i+i> cos t) of degree n. We use inequality (29) with e = — 
and the points 9 3 = —7r, j — 0,..., n — 1 for the first summation, as well as e = ^ and the 
points 9\ = —h-f 77 ", j — 0,..., n — 1 for the second summation. We get 


£ W A \P(A)\’ 


'n+p —1 
2 


< 


X / ( 2 - h*,o)|^(4 +P > c °s^)| r dt 


4lGLD„ 


i=0 


(r + l)e (1 + jh) 

2(n + p)7r 

n+p—1 

+ X J ( 2 — d2i+l,n+p)\P(z2^ P i,COSt)\ r dt J . 


i=0 


Now, on the right hand side, we change the sums with the integrals and apply inequality 
(29) again for the polynomials P(cos s, cost). This time, the variable t is fixed and we use 
e = with the points 9^ = and 9\ = ^d-7r, i — 0,..., n + p — 1. In this way, we 


2t r 
n+p 

obtain 


XI w A\P(A)\ r < (r + l) 2 e 2 

vAGLD n ,p 

= (r + l) 2 e 2 


47t) ^tt 2 


I 1 V 1 

47T / 7T 2 


/ 0 JO 

1 pi 


|P(coss, cost)| r dsdt 
\P(x,y)\ r 


-1 J- 1 \/l - X 2 y/l - + 


dxdy. 


=IIPIi 


□ 
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As a second auxiliary result, we show the inverse quadrature sum estimate: 


Lemma 9. For 1 < r < oo, the inequality 
1 rl rl 


7T 


1 I \P(x,y)\ r ^=^ r^- — - dxdy<D r J 2 w a \P(A)\ 
-i J-' v' 1-x 2 y/1 -y 2 A ^ n 


( 30 ) 


holds for all P G II^ with a constant D r independent of n and p. 


PROOF. We use a duality argument as described in [3T] (and more generally in |22] ) to 
show the inverse inequality. For 1 < r < oo and the dual parameter r' := -, we have the 

representation 

||P|| r = sup (P,0>. 
seW:|| s || r /=i 

Further, if we introduce the partial sums 


S n , P g(x,y) := ^2 Cijfi(x)tj{y), 

(i,j)e r n , p 

Sn, P g{x,y) : = c i,A( x )Tj(y) c i,j = (g,Ti(x)fj(y)), 

(m)erp p 


(31) 


we 


obtain for polynomials P G 11^: 


l|P||r= SUP ( p ,Sn lP g) 

9&L r ' :\\g\\ r ,=l 

< sup (P A n;P r/) + sup c 0 ,„(P,r„( 2 /)). 

9&L r ' :\\g\\ r ,=l g£L r ':\\g\\ r , = l 

Now, we use the fact that C UiP P = P holds for all polynomials P G 11^ and the represen¬ 
tation (25) of C n p P in the Chebyshev basis. In this way, we obtain by the orthonormality 
of the basis Ti(x)Tj(y): 


|P||r < sup ^2 w a P(A) Sn iP g(A ) + 


sup 


c 0 ,r 


H’-'- 1 .AeLD „, p 


;=i 2 


Y, w A P(A)T n (y A ). 


.4eLD„ 


Applying Hblders inequality for the dual pair r, r' to both sums on the right hand side, we 
get 


|P||r < I £«u|iW sup ( *22 w A \Sn )P g(A)\ r ' 

\ 4 ^ T / II17 IF/ — 1 \ yt^TFl 


1 A&'LDr, 


, -4eLD„ 


Y™APW\' 


k AeLDn 


Q) ,n 

sup — 
l.<?L/=i 2 


E ^A\T n (y A )\ r ' 


,.AeLD„ 
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Now, for the first term in the above inequality, we use the forward quadrature sum estimate 
proven in Lemma 8} For the second term, we use the fact that T n (y) < v^2 and, in a second 
step, again the Holder inequality for (g,T n (y)). Then, we obtain 


Mir < 


< 


£ va\p(a)\’ 


. yteLD„ 


£ «uim)r 


. AeLD n 



sup ll^n^llr' + sup 
llsll r /=l IMIr ' =1 


(g,T n (y)) \ 

V2 J 


C;/ ||S n , : 


<p\\ L r ' ^L r ' 



Therefore, if we know that the partial sum operators S UjP are uniformly bounded in the L r 
norm, the proof is finished. To see this, we write S n _ p g as trigonometric partial sum 


S n , p g(cosa,cos (3) = a t j 

r2:=((ij)6Z ! : -fL + M <1 |. (32) 

J i n + p n J 

The index set T^ n ® is the intersection of Z 2 with a rhombus consisting of four congruent right 
triangles with side lengths n-\-p and n. Now, we can adopt a classical result of Fefferman pT6| . 
It states that for 1 < r' < oo the trigonometric partial sum over the rhombic index set T^ 1 ® 
is uniformly bounded in the IF-norm on the 2-torus, i.e. 


II *S'n,p5 , || r' 



a itj e* ai+ M 

(ij)e 



<^||^|| r ,, 


holds for all n,p G N and the constant Cy does not depend on n and p. This proves the 

1 

statement of the Lemma with D r = + l) r . □ 


Remark 6. The combination of forward 


and inverse quadrature sum estimate (30) is 


also referred to as Marcinkiewicz-Zygmund inequality, see 1221 EU and the references therein. 
The idea for the proofs of Lemma [8] and Lemma [9] is taken from [31] where similar results 


are shown for the Xu points. The detailed elaboration of the proofs of the inequalities (28) 
and (30) was necessary in order to guarantee that the constants in (28) and (30) do not 
depend on the parameter p. 


Theorem 10. Let 1 < r < oo, / e C ([—1, l] 2 ) and (p n ) ne ^ a sequence of natural numbers 
such that n and n + p n are relatively prime for all n G N. Then, the Lagrange interpolant 
/~- n ,Pnf converges for n —>■ oo to the function f in the L r -norm, i.e. 

lim \\C n ,p n f ~ f\\r = 0. 
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PROOF. Following the argumentation scheme described in j22], we adopt the inverse quadra¬ 
ture sum estimate of Lemma [ 9 ] to the polynomial C n)Pri f G II^T. In this way, we get for 
1 < r < 00 : 

ii^wii; < d t ;r u>A\f(A)\ r < a-ii/il. 

AgT l 

• AVC - L n,p n 

Now, for an arbitrary polynomial P G II 2 , we have the identity C n>Pn P = P and therefore 

II Cn, P J - f\\r < || Cn, Pn (f ~ P) Ik + 11^ " /||r < (1 + D r )\\P - /11 oo - 

Since polynomials are dense in C([— 1, l] 2 ), we immediately get mean convergence of the 
Lagrange interpolant for 1 < r < oo. The convergence for r = 1 follows analogously using 
the estimate 

IKWIli < I|£„, p „/I|2 < AII/IU 

□ 


Finally, we consider the Lebesgue constant related to the interpolation problem (|19j). It is 
given as the operator norm of the interpolation operator C n , p in the space (C([—1,1] 2 ), || ■ ||oo): 


Ap. 7, . 




max > \L A (x,y)\. 

(*,i/)e[-i,i ] 2 ' 


.AeLD. 


For p — 1, it is known that A njl grows as 0(ln 2 n) |3, IS]. The next theorem states that 
similar behavior is true for general p. 

Theorem 11. The Lebesgue constant A ntP is bounded by 

D a ln 2 (n) < A n>p < C A ln 2 (n + p). 

The constants C\ and D\ do not depend on n and p. 


PROOF. Using formula (24), we get for the Lebesgue constant A. 


n,p • 


A H)P = max 




£ > Wa 

1,1] 2 


.AeLD 


K n , p ( x ,y,XA,yA) + g T n{y) T n{yA) 


< max V w A \K n Ax,y,x A ,y A )\ + l. 


Now, by Lemma [8j we get 

A ntP < Ci max \\K n . p (x,y, -, -)||i + 1. 
U,y)e[-i,i ] 2 
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With the coordinate transforms x = cos a and y = cos/3 we transfer the above norm in a 
trigonometric setting on the 2-torus and obtain 


A <r 

A n, P < max — 

0<a,/3<27r 47 T Z 


[0,27r) 2 


e K ia '+jP') cos(ia) cos(j/3) 


da'dfl' + 1 


1 47T 2 


[0,27t) 2 






dc/d/3' + 1, 


where Tj^ 1 ® denotes the rhombic index set defined in (32). The integral on the right hand side 


corresponds to the L 1 -norm of the Dirichlet kernel with respect to the rhombic summation 
area rj/ 1 ®. By a result of [32] . this norm is bounded by C In (2n + 2 p) with a constant C 
independent of n and p. Thus, we get 

A n p < C\C ln 2 (2n + 2 p) + 1 < (4 C\C + 1) In 2 [n + p) 


and the upper estimate is proven for C\ = ACiC + 1. 

For the lower estimate of the Lebesgue constant, we proceed similar as in a proof given 
for the Padua points [[13]. Since C n . p f is a linear projection of C([—1, l] 2 ) onto If 2 ;^, we get 
by [ 251 Theorem 2.3] the following estimate: 

A n ,p= Slip [IAhp/Hoo > \ Slip ||5^p/Hoo, 

ll/ll=o<l 2 ll/l|oo<l 


where f = (f,K^ p (x,y;-,-)) denotes the partial sum operator given in (J3l]). We note 
that in [29] this inequality is only proven for projections onto fl 2 = If 2 '\. However, a 
straightforward modification of the proof in [29] yields the respective result for general p. 
As a consequence of the Riesz representation theorem (cf. [141 IV.6.3]) we further obtain 


1 

A n,p Pa 2 


max 


\ K n,p( x ,y'r ,-)Hi > X |pW iP (l, 1; 


Using the coordinate transform x = coscc, y = cos/3, we transfer the above norm in the 
trigonometric setting 



dad/3 


with the index set T^’p ng = T^ 1 ® U (0, ±n). Thus, on the right hand side we have again the 
L 1 -norm of a Dirichlet kernel based on a rhombic summation index rA* ng . Since a sphere 
with radius ^ fits into the rhombus with diagonal lengths 2 n and 2n + p we can adopt a 
result of [33] . It states that the Id-norm of the Dirichlet under investigation is then bounded 
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from below by D ln 2 (/h + 1) with a constant D independent of n and p. Therefore, we get 
for Da = f: 

A»,p > y In 2 +lj>j ln 2 (n) = Da ln 2 (n). 

□ 


For a function / G C([— 1, l] 2 ), the best approximation E np (f ) of / in the polynomial 


space II^’p is given 


as 


E n , P (f) : = mi?, 1 , ll/--P||oo- 


ReiF 


For u, v > 0, the modulus of continuity cu(f;u,v) is defined as (cf. PUJ section 3.4.1]) 
<*>(/;«,u):= SU P SU P l/(^i,j/i) - f(x 2 ,y 2 )\- 

\x 1 -x 2 \<u \yi-V 2 \<v 

xi,®2e[-l,l] 2/1,2/2 £[—1,1] 

Using these standard tools from constructive approximation theory together with the esti¬ 
mate of the Lebesgue constant in Theorem 0 we obtain the following error estimates. 


Corollary 12. For any continuous function f G C([— 1, l] 2 ), we have 

\\f ~ £ n , P f Woo < (C\ ln 2 (n + p) + 1 )E nj ,(f). 

If |^, G C([— 1, l] 2 ) for given r,s G N 0 , we further have the estimate 


(33) 


11/ - A^/lloo < C In 2 (n + p) 


U ( dS ’ n+p ’ ^ ) W r>. 


(n + p) r 


+ 


ail-o I 


rr 


(34) 


with a constant C independent of n and p. 

PROOF. We denote by P* the best approximating polynomial of / in II 2 ’/, i.e. \\P* — /||oo — 


E n ,p(f )■ Since C np P* = P*, the estimate of the Lebesgue constant in Theorem 11 leads to 
the following estimate: 

||/ - Wlloo < 11/ - P* lloc + II Cn,pP* - £„,p/||oo 

< (A n ,p + 1)11/ - Ploo = (C A ln 2 (n + p) + l)Pn,p(/). 

Since the polynomial space IT/ = span {Tj(x)Tj(y) : 0 < i < 24^,0 < j < |} is contained 
in n 2 ;/ , we have E np (f) < \\f — Q*\\oc where Q* denotes the best approximating polynomial 
in IT/. Now, a multivariate version of Jackson’s inequality (cf. [5U| section 5.3.2]) gives 


3 r f. 1 


,o 




II f _ /0*|| I 2 ^ydx 7 ' ’n+p’®) , 2 n J \ / P<omax{r,s}+l / ^ V ’ ^+P ’*// . ,W ’ ™ 

11/ W lloo S O I (n+p) r + n s l ( n+p) r + 


In the second inequality, we used the semi-additivity of the modulus c a(f;h x ,h y ). This 
inequality together with (33) yields (34). □ 
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Remark 7. If the function / 6 C([—l, l] 2 ) satisfies the Dini-Lipschitz-type condition 


lim In z {n + p n ) u ( /; —-—, - ) = 0, 
n-> OO \ n + p n n 


inequality (34) guarantees the uniform convergence 


lim 

n—>• oo 


£n,p n f - 


= 0. 


For the Padua points, the result of Corollary [l2| can be found in m The upper estimate for 
the Lebesgue constant of the Padua points proven in [3] is more accurate compared to the 
estimate in Theorem [lTj Further, more recent versions of Jackson’s inequality in a general 


multidimensional setting can be found in [T] and the references therein. 

Remark 8. With slight adaptations of the respective proofs, the convergence results in 
Theorem 10 and Corollary [12] as well as the estimate of the Lebesgue constant in Theorem 
11 can be shown also for the interpolation schemes of the non-degenerate Lissajous curves 
considered in HE In particular, these results confirm the numerical tests given in 


6. Numerical experiments 


Finally, we illustrate numerically the effect of different values of the parameter^, n e N, 
on the Lebesgue constant and the convergence behavior of the polynomial interpolation 
schemes. In particular, we will see that larger values of p n can have advantages when 
approximating functions on anisotropic domains or functions with anisotropic smoothness. 

We investigate first the Lebesgue constant A niPri for the three different parameters p n € 
{1, n + 1, [^n\n + 1}. The first choice p n = 1 leads to respective results of the Padua points 
and can be compared to the numerical experiments given in P!l5j. In Figure |4j(a) the values 
A„. jPn are illustrated for 1 < n < 50. For a better comparison of the values, also the functions 
fi(n) = ln(n) 2 /2 + 2 and / 2 (ri) = In(ny^n) 2 /2 + 4 are plotted in Figure [l|a) as a lower and 
an upper benchmark, respectively. In Figure |4j(b) the Lebesgue constants are plotted with 
respect to the number |LD„ iPn | of interpolation points. The logarithmic growth of A niP „ as 
estimated in Theorem [TT] is clearly visible in Figure [4j Further, the numerical experiments 
indicate a slight growth of A njM with respect to an increasing parameter p n . The best results 
are obtained for p n — 1, i.e. for the Padua points. 

To evaluate the convergence of the interpolation polynomials C n , Pn f to a continuous 
function /, we use the three test functions 


fi(x,y) 

h(x,y ) 
h (x,y) 


(5-10x) 2 (5-10y) 2 (5-10x) 2 +(5-10y) 2 

— e 2 + 0.75e 2 + 0.75e 2 

. . (y-0.5) 2 

= (1[0.25,0.75](t)(T — 0.25) + ll[0.75,oo)(^)j 6 2 , 

= f2(y,x), 


where 11 j(x) denotes the indicator function of an interval J C M, i.e. tj(x) — 1 if x e J 
and llj(a;) = 0 otherwise. The function f\ is taken from the test set in (271 and is smooth, 
whereas / 2 and / 3 have two discontinuities in the partial derivatives ^ and respectively. 
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(a) A n>Pri for n = 1,... 50. (b) A n>Pri with respect to the number |LD„ iPii |. 


Figure 4: The Lebesgue constant A„ jPn for the parameters p n £ {l,n + 1, |_v7ijn + !}• 


As described in M. we use affine mappings of the square [—1,1] 2 to calculate the interpo¬ 
lation polynomials £ n , P „f on the two rectangular domains Q] = [0, l] 2 and fi 2 = [0, 2] x [0,1]. 
As parameters p n we consider again the cases p n G {1, n + 1, \_y/n\n + 1}. The maximal er¬ 
ror between / and C n , Pn f is computed on a uniform grid of 100 x 100 and 200 x 100 points 
defined in fR and 0 2 , respectively. 

The approximation errors for the function fi on the domains fli and 0 2 are displayed in 
Figure [5|a) and (o|b), respectively. While all considered interpolation polynomials converge 
to fi as n —» oo, the decay of the approximation error on the two domains depends on the 
choice of the parameter p n . On the square fR the best results are obtained for the parameter 
p n — 1, whereas the parameter p n — n + 1 seems to produce better adapted interpolation 
polynomials for the anisotropic domain fi 2 . 

Finally, we investigate the influence of p n in the case that the given function is smooth 
in only one of the variables. As test functions, we consider / 2 and / 3 on Op. Since both 
are Lipschitz on fR, the criterion in Remark 7 is satisfied for all three choices of p n . The 
convergence of the respective interpolation schemes is visible in Figure [6j We observe that 
for / 2 the choice p n = |_\/ffil n + 1 produces the best result, whereas for the function / 3 the 
choice p n = 1 leads to a faster decay with respect to the number of interpolation points. 

This numerical result reflects the theoretical estimate given in (34) which contains an 
interplay between two directional moduli of continuity. Since / 2 is smooth with respect to 
y and nonsmooth with respect to x, the modulus o;(/ 2 ; 0, is typically smaller than the 
modulus o;(/ 2 ; n + p , 0) if p n is small. However, using large values for the parameter p n , this 
disparity between the two directional moduli can be compensated and a better approxima¬ 
tion quality can be achieved. For the function / 3 the roles of x and y are interchanged. In 
this case, larger values of p n have no essential effect on the approximation quality. 
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Number of interpolation points 
(cl) ||^3n,p n ^l /l ||oo Oil Oi. 



(b) ||£ n , p „/i - /illoo on n 2 . 


Figure 5: The approximation error \\C n ^ Prl fi — /i||oo ° n the domains fii = [0, l] 2 and tt 2 = [0, 2] x [0,1]. 
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